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DETAILED  MODELLING  OF  COMBUSTION: 

A  NON-INTERFERING  DIAGNOSTIC  TOOL 

I.  Introduction 

Detailed  modelling,  or  numerical  simulation,  provides  a  method  we 
can  use  to  study  complex  reactive  flow  systems.  [1_]  Using  this  tech¬ 
nique,  predictions  about  the  behavior  of  a  physical  system  are  obtained 
by  solving  numerically  the  multi-fluid  conservation  equations  for  mass, 
momentum,  and  energy.  Since  the  success  of  detailed  modelling  is 
coupled  to  the  ability  to  handle  an  abundance  of  theoretical  and  numeri¬ 
cal  detail,  this  field  has  matured  in  parallel  with  the  Increase  in 
size  and  speed  of  computers  and  sophistication  of  numerical  techniques. 

It  is  important  to  distinguish  between  empirical,  phenomenological 
and  detailed  models.  Empirical  models  are  constructed  from  data  ob¬ 
tained  by  experiments,  summarized  in  analytical  or  numerical  form,  and 
subsequently  tested  against  proven  theoretical  laws  or  other  data. 
Phenomenological  models  are  extrapolations  from  theory  based  on  our 
physical  intuition  which  must  be  tested  against  experimental  data.  The 
shortcomings  of  the  empirical  models  lie  in  their  limited  range  of 
validity,  while  phenomenological  models  become  more  tenuous  as  they 
approach  the  complexities  of  real  physical  systems. 

Detailed  models  usually  contain  parts  which  may  be  empirical  or 

phenomenological  in  origin.  However,  detailed  modelling  attempts  to 

overcome  these  shortcomings  by  incorporating  theoretical  detail  rich 

enough  to  approximate  reality;  detail  far  richer  than  could  be  summarized 

in  any  succinct  analytical  model,  yet  more  theoretically  sound  chan 

standard  phenomenological  models  or  empirical  fits. 

Manuscript  submitted  November  6,  1980. 
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The  purpose  of  this  paper  is  to  present  an  overview  and  an 
assessment  of  the  current  state  of  detailed  modelling  as  applied  to 


combustion  systems.  We  will  attempt  to  familiarize  the  reader  with  the 
goals,  terminology  and  inherent  problems  in  modelling  combustion  systems. 

The  emphasis  is  not  on  presenting  a  full  set  of  complicated  multi-fluid 
equations  or  on  explaining  the  numerical  algorithms  required  to  solve 
the  governing  equations.  Instead  we  hope  to  impart  a  sense  of  the 
power  and  role  of  detailed  modelling,  an  understanding  of  why  physical 
insight  must  be  built  into  numerical  algorithms , and  an  indication  of 
how  to  test  these  models  at  every  stage  of  construction  against 
both  theory  and  experiment. 

Figure  1  depicts  the  role  of  both  analytical  and  numerical  modelling 
with  respect  to  experiments.  Increased  accuracy  and  reliability  are  re¬ 
quired  in  proceeding  from  the  evaluation  of  concepts  to  engineering  design. 
The  "calibration"  of  our  understanding  has  been  given  the  pivotal  central 
location.  From  experimental  observation  and  approximate  theoretical  models 
we  can  postulate  quantitative  physical  laws  which  we  expect  an  effect  to 
obey.  These  "laws"  can  be  tested  against  reality  by  incorporating  them  in 
a  detailed  model  which  makes  quantitative  predictions  for  series  of 
experimental  measurements. 

A  computer  simulation  using  a  detailed  model  is  similar  to  an 
experiment  in  that  it  will  not  give  simple  functional  forms  among 
physical  variables.  Each  calculation  is  like  a  unique  experiment  per¬ 
formed  with  one  set  from  an  infinity  of  possible  sets  of  geometric, 
boundary,  and  initial  conditions.  Just  as  valid  results  can  be 
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extracted  from  an  experiment  only  through  an  understanding  of  the 
effects  and  limitations  of  the  instruments  used  in  collecting  the  data, 
results  obtained  using  detailed  modelling  must  be  examined  in  the  light 
of  the  limitations  inherent  in  its  tools,  both  analytical  and  numerical. 
The  first  section  of  this  paper  will  therefore  deal  with  an  exposition 
of  the  problems  inherent  in  detailed  modelling  of  combustion  systems 
so  that  as  we  proceed  we  have  a  healthy  respect  both  for  the  magnitude 
of  the  problems  and  the  limitations  of  our  methods. 

The  next  section  will  concentrate  on  the  choice  of  numerical  algo¬ 
rithms  used  in  the  models.  This  process  corresponds  to  the  construction 
and  design  of  experimental  apparatus  which  must  reflect  a  good  knowledge 
of  the  physics  the  experiment  is  to  study.  Modelling  combustion  systems 
has  its  own  particular  problems  because  of  the  strong  interaction  be¬ 
tween  the  energy  released  from  chemical  reactions  and  the  dynamics  of 
the  fluid  motion.  Release  of  chemical  energy  generates  gradients  in 
temperature,  pressure,  and  density.  These  gradients,  in  turn,  influence 
the  transport  of  mass,  momentum,  and  energy  in  the  system.  On  a  large 
scale,  the  gradients  may  generate  vorticity  or  effect  the  diffusion  of 
mass  and  energy.  On  a  more  microscopic  scale,  they  are  the  origin  of  the 
turbulence  which  drastically  affects  macroscopic  mixing  and  burning 
velocities.  In  modelling  shocks,  detonations,  or  flame  propagation, 
time  and  space  scales  of  interest  can  span  as  many  as  ten  orders  of 
magnitude.  Thus  to  obtain  adequate  resolution,  the  numerical  methods 
must  be  computationally  fast  as  well  as  accurate.  Methods  must  be 
develooed  which  rely  on  asymptotic  solution  techniques  to  follow  short 
time  and  space  scale  phenomena  on  a  macroscale.  It  is  in  this  aspect 


that  detailed  modelling  most  closely  approximates  experiment.  If  our 
numerical  apparatus  cannot  resolve  the  basic  controlling  physical  pro¬ 
cesses,  no  meaningful  calculations  can  be  made  of  their  effects. 

Although  detailed  modelling  does  not  directly  provide  the  types 
of  useful  analytic  relationships  which  guide  our  intuition  and  allow 
us  to  make  quick  estimates,  it  gives  us  the  flexibility  to  evaluate  the 
importance  of  a  physical  effect  by  simply  turning  it  off  or  on  or 
changing  its  strength.  The  model  can  also  be  used  to  test  the  sensiti¬ 
vity  of  the  physical  results  to  independent  theoretical  approximations. 
Those  analytic  results  which  are  available  are  valuable  in  bench¬ 
marking  the  model  in  various  limits.  A  series  of  tests  which  compare 
analytic  results  to  numerical  simulations  may  calibrate  the  simulation 
before  it  is  compared  to  experiments  or  used  for  extrapolation.  Con¬ 
versely,  a  well-tested  detailed  model  serves  as  a  very  useful  means 
of  calibrating  unknown  parameters  and  form  factors  in  approximate 
theories. 

The  last  two  sections  of  this  paper  will  discuss  this  interplay 
between  detailed  modelling  and  theory  and  experiment.  The  third  section 
describes  how  a  model  must  be  tested  in  various  limits  for  physical 
consistency  to  insure  its  accuracy.  The  specific  example  chosen  here 
is  a  comparison  between  an  analytic  solution  and  a  detailed  numerical 
simulation  of  a  premixed  laminar  flame.  The  last  section  shows  how  a 
comparison  between  models  results  and  experiments  can  be  used  to 
calibrate  the  model  and  to  guide  further  experiments.  The  example 
chosen  is  a  calculation  of  flow  over  an  immersed  object  which  is  com¬ 
pared  to  both  experimental  and  theoretical  results. 
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II.  Problems  in  Modelling  Reactive  Flows 


Errors  and  confusion  in  modelling  arise  because  the  complex  set  of 
coupled,  nonlinear,  partial  differential  equations  are  not  usually  an 
exact  representation  of  the  physical  system.  As  examples,  first  con¬ 
sider  the  input  parameters,  such  as  chemical  rates  or  diffusion  coeffi¬ 
cients.  These  input  quantities,  used  as  submodels  in  the  detailed  model, 
must  be  derived  from  more  fundamental  theories,  models  or  experiments. 
They  are  usually  not  krown  to  any  appreciable  accuracy  and  often  their 
values  are  simply  guesses.  Or  consider  the  geometry  used  in  a  calcula¬ 
tion.  It  is  often  one  or  two  dimensions  less  than  needed  to  completely 
describe  the  real  system.  Multidimensional  effects  which  may  be  impor¬ 
tant  are  either  crudely  approximated  or  ignored.  This  lack  of  exact 
correspondence  between  the  model  adopted  and  the  actual  physical  system 
constitutes  the  basic  problem  of  detailed  modelling.  This  problem, 
which  must  be  overcome  in  order  to  accurately  model  transient  combustion 
systems,  can  be  analyzed  in  terms  of  the 

•  multiple  time  scales, 

•  multiple  space  scales, 

•  geometric  complexity,  and/or 

•  physical  complexity 
of  the  systems  to  be  modelled. 

The  first  class  of  problems  arises  as  the  result  of  trying  to 
represent  phenomena  characterized  by  very  different  time  scales.  In 
ordinary  flame  and  detonation  problems  these  scales  range  over  many 
orders  of  magnitude.  When  phenomena  are  modelled  that  have  character¬ 
istic  times  of  variation  shorter  than  the  timestep  one  can  afford,  these 
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phenomena  are  usually  called  "stiff."  Sound  waves  are  stiff  with  re¬ 
spect  to  the  timestep  one  wishes  to  employ  wh,en  modelling  a  subsonic 
flame  speed.  Many  chemical  reaction  rates  are  stiff  with  respect  to 
convection,  diffusion,  or  even  sound  wave  timestep  criteria.  Two 
rather  distinct  modelling  approaches,  global  implicit  and  timestep- split 
asymptotic,  have  been  developed  to  treat  these  temporally  stiff 
phenomena.  These  two  approaches  are  briefly  described  later  in  this 
paper. 

The  second  class  of  problems  involves  the  huge  disparity  in  space 
scales  occurring  in  combustion  problems.  To  model  the  steep  gradients 
at  a  flame  front,  a  cell  spacing  of  10“ 3  cm  or  smaller  might  be  re¬ 
quired.  To  model  convection,  grid  spacings  of  1  to  10  cm  might  be 
adequate.  Complex  phenomena  such  as  turbulence  which  occur  on  inter¬ 
mediate  spatial  scales  present  a  particular  modelling  problem.  It 
would  be  a  pipedream  to  expect  a  numerical  calculation  to  faithfully 
reproduce  physical  phenomena  with  scale  lengths  shorter  than  a  cell 
size.  Therefore,  to  calculate  realistic  profiles  of  physical  variables, 
a  certain  cell  spacing  is  required  to  obtain  a  given  accuracy.  Choos¬ 
ing  a  method  which  maximizes  accuracy  with  a  minimum  number  of  grid 
points  is  a  major  concern  in  detailed  modelling. 

The  third  set  of  obstacles  arises  because  of  the  geometric  com¬ 
plexity  associated  with  real  systems.  Most  of  the  detailed  models 
developed  to  date  have  been  one-dimensional,  but  this  gives  a  very 
limited  picture  of  how  the  energy  release  affects  the  hydrodynamics. 

Even  though  many  processes  in  a  combustion  system  can  be  modelled  in 
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one-dimension,  there  are  others,  such  as  boundary  layer  growth,  or  the 
formation  of  vortices  and  separating  flows,  which  clearly  require  at 
least  two-dimensional  hydrodynamics.  Real  combustion  systems  are  at 
least  two-dimensional,  with  unusual  boundary  conditions  and  internal 
sources  and  sinks.  However,  even  with  sixth  generation  parallel  pro¬ 
cessing  computers  available,  what  can  be  achieved  with  two-dimensional 
detailed  models  is  still  limited  by  computer  time  and  storage  require¬ 
ments  . 

In  the  current  state-of-the-art,  one-dimensional  models  can  best 
be  used  to  look  in  detail  at  the  coupling  of  a  very  large  number  of 
species  interactions  in  a  geometry  that  is  an  approximation  to  reality. 
Processes  such  as  radiation  transport,  turbulence,  or  the  effects  of 
heterogeneity  of  materials  can  be  included  either  as  empirically  or 
theoretically  derived  submodels.  Two-  and  three-dimensional  models  are 
best  used  to  study  either  gross  flow  properties  or  detailed  radiation 
transport.  In  these  latter  models,  the  chemical  reaction  scheme  is 
usually  quite  idealized  or  parameterized. 

The  final  set  of  obstacles  to  detailed  modelling  concerns  physical 
complexity.  Combustion  systems  usually  have  many  interacting  species. 
This  leads  to  sets  of  many  coupled  equations  which  must  be  solved 
simultaneously.  Complicated  ordinary  differential  equations  describing 
che  chemical  reactions  or  large  matrices  describing  the  molecular  diffu¬ 
sion  process  are  costly  and  increase  calculation  time  orders  of  magni¬ 
tude  over  idealized  or  empirical  models.  Table  1  lists  some  of  the 
major  chemical  and  physical  processes  which  have  to  be  considered  for 
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an  accurate  description  of  a  complicated  combustion  system.  Multi¬ 
phase  processes  such  as  surface  catalysis  and  soot  formation  can  be 
important  even  when  we  are  primarily  interested  in  gas  phase  combustion. 
For  most  interesting  systems,  one  finds  that  the  basic  chemical  reaction 
scheme,  the  individual  chemical  rates,  the  optical  opacities,  or  the 
effects  of  surface  reactions  are  not  well  known.  Before  a  model  of  a 
whole  combustion  system  can  be  assembled,  each  individual  process  must 
be  separately  understood  and  modelled.  These  submodels  are  either 
incorporated  into  the  larger  detailed  model  directly  or,  if  the  time 
and  space  scales  are  too  disparate,  they  must  be  fit  in  phenomenologi¬ 
cally.  For  example,  diffusion  and  thermal  conductivity  between  a  wall 
and  the  reacting  gas  can  be  studied  separately  and  then  incorporated 
directly  into  a  detailed  combustion  model.  Turbulence,  however,  can  be 
modelled  on  its  own  space  scales  only  in  idealized  cases.  These  more 
fundamental  models  must  be  used  to  develop  phenomenological  models  for 
use  in  the  macroscopic  detailed  model.  Resolution  and  computational 
cost  prevent  incorporating  the  detailed  turbulence  model  directly. 

Table  1 

Fundamental  Processes  in  Combustion 

gas  phase  multi-phase 

Chemical  kinetics 
Hydrodynamics- laminar 
Thermal  conductivity,  viscosity 
Molecular  diffusion 
Thermochemistry 
Hydrodynamics- turbulent 
Radiation 

Nucleation  t 

Surface  Effects  ! 

Phase  Transitions  *  * 

(Evaporation,  condensation. . . ) 
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Often  there  are  cases  where  the  submodels  are  poorly  known  or 
misunderstood,  and  this  may  cause  the  modellers  no  end  of  grief.  A 
typical  example  is  shown  in  Fig.  2  which  was  provided  by  David  Garvin 
at  the  U.  S.  National  Bureau  of  Standards.  The  figure  shows  the  rate 
at  300°K  for  the  reaction  HO  +  0^  HO  +  0^  as  a  function  of  the  year 
of  the  measurement.  Me  note  with  amusement  and  chagrin  that  if  we  were 
modelling  a  kinetics  scheme  which  incorporated  this  reaction  before 
1970,  the  rate  would  be  uncertain  by  five  orders  of  magnitude!  Similar 
tales  of  horror  also  exist  for  thermochemical  data.  In  combustion  and 
atmospheric  physics,  detailed  models  have  been  used  to  bound  and  test 
the  importance  of  chemical  rates.  Some  of  these  combustion  applica¬ 
tions  will  be  described  below. 

In  order  to  illustrate  how  the  problems  caused  by  the  requirements 
of  temporal  and  spatial  resolution  and  geometric  and  physical  complexity 
are  translated  into  computational  time,  we  have  chosen  to  analyze  a 
gedanken  flame  experiment.  Consider  a  closed  tube  one  meter  long  which 
contains  a  combustible  gas  mixture.  We  wish  to  calculate  how  the 
physical  properties  such  as  temperature,  species  densities,  and  posi¬ 
tion  of  the  flame  front  change  after  the  mixture  is  ignited  at  one  end. 
The  burning  gas  can  be  described,  we  assume,  by  a  chemical  kinetics 
reaction  rate  scheme  which  involves  some  tens  of  species  and  hundreds 
of  chemical  rates,  some  of  which  are  "stiff".  We  will  assume  one¬ 
dimensional  propagation  along  the  tube.  Boundary  layer  formation  and 
turbulence  will  be  ignored.  We  further  assume  that  the  flame  front 
moves  at  an  average  velocity  of  100  cm/sec. 


10 


A  MODELLER’S  NIGHTMARE 
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YEAR 

Fig.  2  -  Measured  values  of  the  rate  constant  for  HO  +  03  HO2  +  O2  as  a  function  of  the  yea 
measurement.  The  arrows  with  overbars  and  underbars  indicate  measured  upper  and  lower  bounds 
respectively.  (Courtesy  of  David  Garvin,  U.  S.  National  Bureau  of  Standards). 


Table  II  summarizes  the  pertinent  time  and  space  scales  in  this 


problem.  Assuming  the  speed  of  sound  is  105  cm/ sec,  a  timestep  of 
about  10  ^  sec  sould  be  required  to  resolve  the  motion  of  sound  waves 
bouncing  across  the  chamber.  Chemical  timescales, as  mentioned  above, 
are  about  10-6  sec.  This  number  may  be  reduced  drastically  if  the 
reaction  rates  or  density  changes  are  very  fast.  It  takes  a  sound 
wave  about  10~3  seconds  to  cross  the  1  meter  system  and  it  takes  the 
flame  front  about  one  second  to  cross.  We  further  assume  that  the 
flame  zone  is  about  10-2  cm  wide  and  that  it  takes  grid  spacings  of 
10* 3  cm  to  resolve  the  steep  gradients  in  density  and  temperature.  In 
those  portions  of  the  tube  on  either  side  of  the  flame  front,  we  assume 
that  1  cm  spacings  are  adequate. 

Table  II 

Important  Scales  in  Gedanken  Flame  Calculation 
Timescales  Spacescales 


At 

sec 

Ax 

cm 

Sound  Speed 

io- 9 

Flame  Resolution 

10“ 

Chemistry 

10"6 

Flame  Zone 

10“ 

Sound  Transit 

Time 

10- 3 

Diffusion  Scale 

10“ 

Flame  Transit 

Time 

1 

Convective  Scales 

10 

System  Size 

100 

Vf 

-  100 

cm/ sec 

To  estimate  the  computational  expense  of  this  calculation,  we  use 
10~3  seconds  of  computer  time  as  a  reasonable  estimate  of  the  time  it 
takes  to  integrate  each  grid  point  for  one  timestep  (a  single  point- 
step)  .  This  estimate  includes  a  solution  of  all  of  the  chemical  and 
hydrodynamic  equations  and  is  based  on  a  detailed  model  of  a  hydrogen- 
oxygen  flame  problem  optimized  for  a  parallel  processing  computer. 


Figure  3  shows  the  information  in  Table  II  cast  into  a  graph  of 
space  versus  time.  Since  the  scales  are  logarithmic,  a  calculation  of 


the  number  of  point-steps  and  then  of  the  needed  computer  time  requires 
exponentiation.  Thus  it  appears  that  3000  years  of  computer  time  is 
required  to  calculate  the  101**  point-steps  involved  in  representing 
the  finest  resolved  space  and  time  scales! 

Of  course  this  is  unacceptable.  Ideally  such  a  simple  calculation 
should  take  about  100  seconds.  What  are  needed  are  numerical  algorithms 
which  have  the  resolution  in  time  and  space  only  where  it  is  required. 
Furthermore,  these  algorithms  should  be  optimized  to  take  advantage  of 
what  is  known  about  the  physics  and  chemistry  of  the  problem.  This 
will  be  discussed  further  below  where  it  is  shown  how  the  application 
of  various  numerical  algorithms  can  be  used  to  reduce  this  flame 
system  to  a  tractable  computational  problem. 

Turbulence  is  one  of  the  outstanding  problems  of  reactive  flow 
modelling  and  is  another  excellent  example  of  the  difficulty  we  have 
in  resolving  highly  disparate  time  and  space  scales.  Our  understanding 
and  eventual  ability  to  predict  the  complicated  interactions  occurring 
in  turbulent  reactive  flow  problems  is  imperative  for  many  combustion 
modelling  applications.  The  presence  of  turbulence  alters  mixing 

and  reaction  times  and  heat  and  mass  transfer  rates  which  in  turn  modi¬ 
fy  the  local  and  global  dynamic  properties  of  the  system.  What  we 
need  to  resolve  these  problems  are  accurate  yet  compact  phenomenologi¬ 
cal  turbulence  models  which  can  be  used  to  describe  realistic  combustor 
systems,  open  flames,  and  other  turbulent  reactive  flows  confidently 
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Fig.  3  -  Space  and  time  scales  In  the  gedanken  flame  calculation.  A  naive  direct  solution  of  the 
problem  could  take  3000  years  of  computer  time.  The  calculation  should  be  possible  In  100  seconds 


and  efficiently.  These  computational  models  must  decouple  the  subgrid 
turbulence  and  microscopic  instability  mechanisms  from  calculations  of 
the  macroscopic  flow.  Below  we  list  the  important  properties  of  an 
ideal  turbulence  model  [2J . 

1.  Chemistry-Hydrodynamic  Coupling  and  Feedback 

Explicit  energy  feedback  mechanisms  from  mixing  and  reactions  to 
the  turbulent  velocity  field  and  the  macroscopic  flow  must  be  formu¬ 
lated.  The  "laminar"  macroscopic  flow  equations  contain  phenomeno¬ 
logical  terms  which  represent  averages  over  the  macroscopic  dynamics 
to  include  the  effects  of  turbulence.  Examples  of  these  terms  are 
eddy  viscosity  and  diffusivity  coefficients  and  average  chemical  heat 
release  terms  which  appear  as  sources  in  the  macroscopic  flow  equations. 
Besides  providing  these  phenomenological  terms,  the  turbulence  model 
must  use  the  information  provided  by  the  large  scale  flow  dynamics 
self-consistently  to  determine  the  energy  which  drives  the  turbulence. 
The  model  must  be  able  to  follow  reactive  interfaces  in  the  macroscopic 
scale. 

2.  Modelling  Onset  and  other  Transient  Turbulence  Phenomena 

The  model  should  be  able  to  predict  the  onset  of  turbulence  in 
initially  laminar  flow  since  bursts  and  other  highly  transient  pheno¬ 
mena  seem  to  be  the  rule  in  reactive  flow  turbulence.  Gradients  in 
density,  temperature,  and  velocity  fields  in  the  reacting  fluid 
drive  the  macroscopic  fluid  dynamic  instabilities  which  initiate  tur¬ 
bulence.  Thus  these  gradients  from  the  macroscopic  calculation  are 
bound  to  be  key  ingredients  in  determining  the  energy  available  to 
drive  the  turbulence. 
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3.  Complicated  Reactions  and  Flow 

The  ideal  turbulence  model  must  deal  with  multiscale  effects  within 
the  subgrid  model.  If  there  is  a  delay  as  velocity  cascades  to  the 
short  wavelength  end  of  the  spectrum  due  to  chemical  kinetics  or 
buoyancy,  for  example,  the  model  must  be  capable  of  representing  this. 
Otherwise  bursts  and  intermittency  phenomena  cannot  be  calculated. 

4 .  Lagrangian  Framework 

An  ideal  subgrid  model  should  be  constructed  on  a  Lagrangian  hydro¬ 
dynamics  framework  moving  with  the  macroscopic  flow.  This  requirement 
reduces  purely  numerical  diffusion  to  zero  so  that  realistic  turbulence 
and  molecular  mixing  phenomena  will  not  be  masked  by  non-physical  numeri¬ 
cal  smoothing.  This  requirement  also  removes  the  possibility  of  masking 
purely  local  fluctuations  by  truncation  errors  from  the  numerical  repre¬ 
sentation  of  macroscopic  convective  derivatives.  The  time-dependent 
(hyperbolic)  Lagrangian  framework  should  also  generalize  to  three  dimen¬ 
sions  as  well  as  resolve  reactive  interfaces  dynamically. 

5.  Scaling 

Breaking  a  calculation  into  macroscopic  scales  and  subgrid-scales 

is  an  artifice  to  allow  us  to  model  turbulence.  The  important  physics 

occurs  continuously  over  the  whole  spectrum  from  k  ,  the  wavenumber 

corresponding  to  the  system  size, to  k  ,  the  wave  number  corres- 

diss. 

ponding  to  a  mean  free  path  of  a  molecule.  Thus  the  macroscopic  and 

subgrid  scale  spectra  of  any  physical  quantity  must  couple  smoothly  at 

k  the  cell  boundary  wave  number.  If  this  number  were  to  be 

cell 

changed,  as  might  happen  if  numerical  resolution  were  halved  or 
doubled,  the  predictions  of  the  turbulence  model  must  not  change. 
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6.  Efficiency 

Of  course,  the  model  must  be  efficient.  The  number  of  degrees  of 
freedom,  required  to  specify  the  status  of  turbulence  in  each  separately 
resolved  subgrid  region  has  to  be  kept  to  a  minimum  for  the  model  to 
be  generally  useable.  The  real  fluid  has  essentially  an  infinite 
number  of  degrees  of  freedom  to  represent  the  state  of  the  gas  in  each 
small  element.  We  would  like  to  be  able  to  do  the  job  with  a  minimal 
number  of  degrees  of  freedom. 

III.  Choosing  an  Algorithm  Based  on  the  Physics  of  the  Problem 

In  reactive  flow  calculations  we  are  concerned  with  two  flow 
regimes  which  depend  on  the  rate  of  energy  release.  When  energy  is 
released  quickly,  shocks  and  detonations  are  formed.  When  energy  is 
released  slowly,  flames  are  formed.  The  former  requires  that  the 
numerical  algorithm  used  follow  the  changes  of  the  system  on  time- 
scales  determined  by  the  speed  of  sound  in  the  material  (Courant 
condition).  If  we  follow  this  same  timescale  in  the  flame  case,  in 
which  the  physical  timescales  of  interest  are  much  larger,  the  cost  is 
exhorbitant.  The  gedanken  flame  calculation  described  above  cost  so 
much  partly  because  we  postulated  the  use  of  an  explicit  algorithm 
based  on  timesteps  determined  by  the  Courant  condition.  For  flame 
calculations,  then,  the  answer  is  to  use  techniques  in  which  the  energy 
conservation  equation  is  converted  to  a  pressure  equation  which  is 
solved  implicitly. 

The  problem  is  basically  that  of  coupling  into  one  calculation  all 
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of  the  pertinent  physical  and  chemical  processes  characteristic  of  a 
combustion  system.  Two  distinct  approaches  have  evolved.  In  the  first 
of  these,  often  called  "global  implicit"  differencing,  the  complete 
set  of  nonlinear  coupled  equations  describing  the  physical  system  of 
interest  is  cast  into  a  simple  finite-difference  form.  The  spatial 
and  temporal  derivatives  are  discretized  and  the  nonlinear  terms  are 
linearized  locally  about  the  solutions  obtained  numerically  at  the 
previous  timestep.  This  process  is  valid  only  when  the  values  of  the 
physical  variables  change  slowly  over  a  timestep.  A  rigorously  correct 
treatment  of  the  nonlinear  terms  reqyires  iteration  and  large  matrices 
must  be  inverted  at  each  timestep  to  guarantee  stability.  In  one 
spatial  dimension  the  problem  usually  appears  as  a  block  tridiagonal 
matrix  with  M  independent  physical  variables  to  be  specified  at 
grid  points.  Then  an  MN^  by  MNx  matrix  must  be  inverted  at  each 
iteration  of  each  timestep.  The  blocks  on  or  adjacent  to  the  matrix 
diagonal  are  M  x  M  in  size  so  the  overall  matrix  is  quite  sparse. 
Nevertheless,  an  enormous  amount  of  computational  work  goes  into 
advancing  the  solution  even  a  single  timestep.  Multidimensional  pro¬ 
blems,  in  this  approach,  lead  to  matrices  which  areMN  N  by  MN  N  in 

x  y  x  y 

two  dimensions  and  MN  N  N  by  MN  N  N  in  three  dimensions.  In  complex 

xyz  xyz 

kinetics  problems  with  no  spatial  variation,  the  M  independent  vari¬ 
ables  are  the  species  number  densities  and  temperature  in  the  homo¬ 
geneous  volume  of  interest.  The  Gear  method  [_3]  is  an  example  of  this 
global  implicit  approach  for  pure  kinetics  problems. 

The  second  approach  is  a  fractional-step  method  we  call  asymptotic 
timestep-splitting .  It  is  developed  by  consideration  of  the 
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specific  physics  of  Che  problem  being  solved.  Stiffness  in  the  gover¬ 
ning  equations  is  usually  handled  "asymptotically"  rather  than  impli¬ 
citly.  The  individual  terms,  including  those  which  lead  to  the  stiff 
behavior,  are  solved  as  independently  and  accurately  as  possible. 
Examples  of  such  methods  include  the  Selected  Asymptotic  Integration 
Method  [4.,5J  for  kinetics  problems  and  the  asymptotic  slow  flow  algo¬ 
rithm  for  hydrodynamic  problems  where  the  sound  speed  is  so  fast  that 
the  pressure  is  essentially  constant  [£,21  • 

The  tradeoffs  between  these  two  approaches  are  clear.  The  implicit 
approach  puts  maximum  strain  on  the  computer  and  minimal  strain  on  the 
modeller.  For  this  method,  convergence  of  the  computed  solutions  is 
easy  to  test  with  improved  temporal  and  spatial  resolution.  Mon¬ 
convergence  of  any  particular  calculation  may  be  hard  to  spot  since 
severe  numerical  damping  has  been  introduced  to  maintain  numerical 
stability  and  positivity.  This  damping  changes  the  desired  profiles 
quantitatively,  although  quickly  detected  qualitative  errors  are  often 
smoothed  out.  Solutions  may  be  wrong  yet  stable. 

In  contrast,  the  asymptotic  approach  puts  minimal  strain  on  the 
computer  but  demands  much  more  of  the  modeller.  The  convergence  of  the 
computed  solutions  is  usually  easy  to  test  with  respect  to  spatial  and 
temporal  resolution,  but  situations  do  exist  where,  for  example,  re¬ 
ducing  the  timestep  can  make  an  asymptotic  treatment  of  a  "stiff" 
phenomenon  less  accurate  rather  chan  more  accurate.  This  follows 
because  the  disparity  of  timescales  between  fast  and  slow  phenomena  is 
often  exploited  in  the  asymptotic  approach  rather  than  tolerated. 
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Furthermore,  the  non-convergence  of  any  particular  solution  is  often 
easier  to  spot  in  the  asymptotic  approach  because  the  manner  of  degra¬ 
dation  is  usually  catastrophic.  In  kinetics  calculations,  lack  of 
conservation  of  mass  or  atoms  signals  inaccuracy  rather  clearly. 

The  asymptotic  approach  usually  leads  to  more  modular  simulation 
models  than  the  global  implicit  approach.  Hydrodynamics,  transport, 
equation  of  state  calculations ,  and  chemical  kinetics  are  tied  neatly 
into  individual  packages.  What  is  even  more  important,  specialized 
techniques  for  enhancing  accuracy  can  be  incorporated  at  each  stage  and 
for  each  physical  phenomenon  being  modelled  separately.  There  is  no 
need  to  use  simpler  methods  which  are  suitable  for  inclusion  into  a 
single  giant  finite  difference  formula.  Since  each  phenomenon  is  treated 
as  an  independent  package,  the  full  spectrum  of  numerical  tricks  is 
applicable. 

These  packages  are  relatively  easy  to  test  individually  and  can  be 
very  sophisticated.  They  can  also  be  used  directly  in  a  number  of 
totally  different  physical  problems  with  little  or  no  change  and  are 
hence  more  flexible  than  equivalent  portions  of  a  global  implicit 
algorithm.  The  price  for  this  flexibility  is  the  need  to  treat  care¬ 
fully  all  the  couplings  between  the  individual  physical  terms  and 
effects.  Using  the  asymptotic  approach  one  cannot  sit  back  and  turn  a 
massive  mathematical  crank  to  get  an  answer. 

At  this  point  the  pros  and  cons  of  the  two  approaches  seem  to 
roughly  counterbalance;  they  have  been  presented  that  way  purposely. 

This  apparent  equity  extends  to  most  accuracy  criteria  as  well.  If 


an  interesting  timescale  is  not  resolved,  neither  solution  method 
can  give  detailed  profiles  of  phenomena  occurring  on  that  scale. 
Similarily,  to  compute  spatial  gradients  accurately  they  must  be  re¬ 
solved  with  enough  spatial  grid  points  in  either  type  of  calculation. 

The  fact  that  the  asymptotic  approach  demands  more  work,  of  the 
physicist  is  counterbalanced  by  the  work,  that  must  be  done  to  reduce  the 
computational  expense  of  using  the  global  implicit  method.  This 
calculational  expense,  above  all  else,  is  the  factor  which  has  caused 
us  to  employ  asymptotic  rather  than  global  implicit  formulations.  For 
example,  solving  a  chemical  kinetics  scheme  for  M  species  requires 
inverting  a  general  matrix  of  size  Mx  M.  This  involves  approximately 
M3  operations.  In  contrast,  the  selected  asymptotic  approach  to  sol¬ 
ving  the  kinetics  equations  generally  scales  as  M.  It  is  one  goal  of 
detailed  modelling  to  be  able  to  include  the  full  details  of  extremely 
complex  kinetics  systems  coupled  to  time-dependent  fluid  dynamics. 

Since  more  complex  problems  can  be  solved  for  the  same  cost  using 
asymptotics,  we  are  willing  to  invest  the  effort  in  the  physics  modules 
and  their  coupling  in  order  to  be  able  to  expand  our  computational 
abilities. 

Using  the  information  discussed  so  far,  we  can  now  return  to  the 
gedanken  flame  calculation  with  the  idea  of  modifying  our  numerical 
methods  in  order  to  reduce  the  computational  cost.  The  goal  is  to 
calculate  the  propagation  of  a  flame  front  across  a  one-meter  tube 
using  a  one-dimensional  geometry  and  given  a  fixed  detailed  chemical 
reaction  rate  scheme. 
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First,  we  recognize  immediately  that  we  are  interested  in  calcula¬ 


ting  a  flame  front  moving  at  less  than  the  local  sound  speed.  Thus 
either  a  slow  flow  approximation  or  any  method  which  treats  pressure 
implicitly  would  eliminate  the  sound  speed  criterion  on  the  timestep. 

By  using  the  asymptotic  slow  flow  technique  described  below  and  still 
assuming  a  uniform  grid  spacing,  the  number  of  point-steps  is  reduced 
from  1014  to  1011.  Thus  Fig.  4  shows  that  the  time  required  for  the 
calculation  is  reduced  from  3000  to  3  years! 

But  this  is  still  atrocious,  and  we  must  now  face  the  problem  of 
eliminating  unnecessary  grid  points.  Adaptive  gridding  is  currently 
a  frontier  in  reactive  flow  modelling.  As  we  have  mentioned  above, 
there  are  no  excellent  techniques.  The  block  on  the  graph  in  Fig.  5 
shows  the  region  spanned  in  the  gedanken  flame  problem  by  an  adaptively 
gridded  calculation.  Here  100  cells  of  1  cm  length  are  used  and  the 
region  surrounding  the  flame  front  is  finely  gridded  with  100  addi¬ 
tional  cells  of  10-3  cm  length.  The  timestep  is  still  governed  by  the 
smallest  cells,  but  by  now  only  200  cells  are  needed  rather  than  105. 
The  saving,  about  a  factor  of  500,  reduces  the  computational  time  to 
2  x  10s  point-steps,  or  about  two  days. 

Finally,  Fig.  6  summarizes  the  computational  expense  of  performing 
the  flame  propagation  problem  using  the  possible,  but  as  yet  unexploited, 
technique  of  adaptive  intermittent  gridding.  The  idea  here  is  that  a 
finely  gridded  region  is  injected  into  the  calculation  at  intermittent 
timesteps.  This  is  done  often  enough  to  update  the  properties  of  the 
finely-spaced  region  which  are  then  used  as  interior  boundary  condi¬ 
tions  for  the  coarsely-spaced  region.  Now  assume  that  100  ceils  are 
needed  to  resolve  the  flame  zone.  Further,  100  short  timesteps  are 
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Fig.  4  -  Using  the  slow  flow  technique,  which  allows  us  to  follow  timescales  larger  than  those  required 
by  an  explicit  solution  of  the  energy  equation,  reduces  the  required  computational  time  to  three  years. 
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Ats  Ate  SOUND  FLAME 

Fig.  6  -  The  as  yet  undeveloped  injected  adaptive  gridding  method  will  reduce  the  computational  time 

to  1100  seconds 


enough  Co  resolve  changes  in  the  flame  zone  brought  about  by  the 
relatively  slowly  changing  outer  boundary  conditions.  During  the  im¬ 
bedded  calculation,  the  flame  front  moves  only  10  of  the  fine  zones, 
which  is  sufficient  to  determine  flame  speed  and  boundary  conditions 
to  be  used  in  the  coarsely  spaced  calculation.  The  imbedded  calcu¬ 
lation  may  be  done  once  in  each  large  cell.  Thus  a  total  of  100  + 

(100) (10)  *  1100  seconds  of  computational  time  is  required  for  the  large 
scale  simulation. 

This  example  has  illustrated  the  importance  of  using  the  appro¬ 
priate  algorithm  motivated  by  considerations  of  the  actual  problem 
that  must  be  resolved.  It  has  further  illustrated  how  much  may  be 
accomplished  by  developing  the  methods  of  adaptive  gridding.  One 
point  that  has  not  been  mentioned,  however,  is  that  much  of  the  cost 
of  a  detailed  reactive  flow  calculation  is  taken  up  by  the  integration 
of  the  ordinary  differential  equations  describing  the  chemical  kinetics. 
Using  the  latest  asymptotic  techniques  improves  the  picture  painted 
above  by  a  factor  of  two  to  four.  But  further  improvements  in  the 
integration  time  of  ODE's  without  sacrificing  accuracy  is  certainly  an 
area  where  development  is  needed. 

IV.  Testing  the  Model  Against  Theoretical  Results 

Analytical  solutions,  while  often  approximate,  are  extremely  useful 
in  providing  functional  relationships  and  generalizing  trends .  Below 
we  show  that  by  comparing  numerical  and  analytic  results,  we  can  gain 
new  insights  into  the  controlling  physical  processes. 

Ignition  of  a  fuel-oxidizer  mixture  occurs  when  an  external  source 
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of  energy  initiates  interactions  among  the  controlling  convective, 
transport  and  chemical  processes.  Whether  the  process  results  in  defla¬ 
gration,  detonation,  or  is  simply  quenched  depends  on  the  intensity, 
duration,  and  volume  affected  by  an  external  heat  source.  Ignition 
also  will  depend  on  the  initial  ambient  properties  of  the  mixture  which 
determine  the  chemical  induction  time  and  the  heat  release  per  gram. 

Thus  ignition  is  a  complicated  phenomena  and  its  prediction  for  a 
specific  mixture  of  homogeneous,  premixed  gases  depends  strongly  on 
input  parameters  which  are  often  very  poorly  known.  A  convenient, 
inexpensive  way  to  estimate  whether  a  mixture  will  ignite  given  a 
heat  source  intensity,  duration,  and  volume  would  be  a  valuable  labora¬ 
tory  tool  and  a  useful  learning  device. 

A  closed  form  similarity  solution  for  the  nonlinear  time-dependent 
slow-flow  equations  has  been  used  as  the  basis  for  a  simple,  time- 
dependent,  analytic  model  of  localized  ignition  which  requires  minimal 
chemical  and  physical  input  [8J .  As  a  fundamental  part  of  the  model , 
there  are  two  constants  which  must  be  calibrated:  the  radii, or  fraction 
of  the  time-dependent  similarity  solution  radius, at  which  the  thermal 
conductivity  and  induction  parameters  are  evaluated.  This  calibration 
is  achieved  by  comparison  with  the  results  of  a  detailed  time-dependent 
numerical  flame  simulation  model  which  is  a  full  solution  of  the 
multi-fluid  conservation  equations.  The  detailed  model  itself  has  been 
checked  extensively  with  respect  to  its  various  chemical,  diffusive 
transport,  and  hydrodynamic  components. 
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The  basic  similarity  solution  is  derived  from  the  slow  flow 


approximation,  characterized  by  (1)  flow  velocities  which  are 
small  compared  to  the  speed  of  sound,  and  (2)  an  essentially  constant 
pressure  field.  The  energy  and  velocity  equations  may  then  be  written 
as 

dP  -lr2  (  . \  i«2 

5£.  o  *  -  yPV  •  v  +  V  •  yNk  <VT  +  S(t)e  ,  (1) 

dt  —  —  —  B  — 

from  which  we  can  derive  an  algebraic  equation  for  7_  •  v.  Here  P  is 
the  total  pressure,  v  is  the  fluid  velocity,  T  is  the  temperature,  y  is 
the  ratio  of  heat  capacities  Cp/Cv  assumed  here  to  be  a  constant  and  < 
is  a  function  of  the  mixture  thermal  conductivity,  X  , 


s  lzi_ 

"  yNk* 


\n(T)- 


(2) 


The  last  term  on  the  right  hand  side  of  Eq.  (1)  is  the  source  term. 

Proper  choice  of  S(t)  ensures  that  a  given  amount  of  energy,  E  ,  is 

o 

deposited  in  a  certain  volume,  —  Rq,  in  a  time,  The  choice  of  this 

Gaussian  profile  allows  us  to  obtain  a  "closed"  form  similarity  solu¬ 
tion.  If  the  fluid  velocity  v  is  then  expanded  such  that 


v(r,t) « v^(t)r 


(3) 


and  spherical  symmetry  is  assumed,  Eq.  (1)  may  be  solved  analytically 
to  obtain 


T(r, t) 


A(t)e 


-k2(t)r2 


(4) 
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and 


P (r,c) 


-A(t)e 


-k2 (t)r2 


(5) 


where  and  are  the  background  temperature  and  density  far  from  the 
heat  source.  Solutions  for  A(t)  and  k(t)  may  then  be  obtained  by  sol¬ 
ving  two  coupled  differential  equations 


dk 

dt 


(6) 


dA 

dt 


lit i_ 

yp 


6<k2A. 


(7) 


Invoking  energy  conservation  yields  the  expression  for  the  velocity 
coefficient  v^,  and  is  effectively  the  first  calibration  in  the  model: 


v 


1 


_S__  F7  (O)-F'  (A)  2  AF'  (A)-F(A) 

3y?x  F(A)  L  ^  F(A) 


(8) 


where  the  function  F(A)  is  defined  as 


oo  — 

F(A)  H  /  4rrx2[l-e'Ae  ]dx. 


(9) 


The  model  requires  one  further  definition  in  order  to  predict 
Ignition.  A  curve  of  chemical  induction  time  as  a  function  of  tempera¬ 
ture  must  be  included  in  order  to  define  the  induction  parameter, 
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(10) 


I(t)  *  L  r  (T(r,  t'  )) 
c 

Ignition  "occurs"  when  I(t)  -  1  in  this  model,  which  is  an  exact  result 
in  the  limit  of  large  heat  source  and  constant  temperature  near  the 
center  of  the  heated  region.  A  simple  analytic  expression  for  t^(T) 
depending  on  three  constants  has  been  derived  and  can  be  calibrated 
using  as  few  as  three  distinct  values  of  at  different  temperatures. 


The  chemical  reaction  scheme  used  in  the  detailed  model  was  used 
to  generate  a  curve  for  tc(T).  The  values  of  thermal  conductivity 
used  in  the  detailed  model  were  used  to  generate  a  function  k.  Then  a 
series  of  comparisons  were  made,  in  which  the  detailed  model  was  con¬ 
figured  in  spherical  symmetry  with  a  Gaussian  energy  deposition. 

We  show  results  for  several  test  cases.  In  one,  R  »  0.1  cm  and 

o 

tq  *  IxlO-4  sec.  The  simple  model  predicts  that  3.3xl04  ergs  is  the  mini¬ 
mum  ignition  energy  and  these  results  agree  well  with  the  simulation 
(Figs.  7  and  8).  Both  models  predict  ignition  at  essentially  the  same 

time  for  a  range  of  input  energies.  In  the  second  example,  R  ■  0.025  cm 

o 

and  t  ■  lxlO”4  sec.  The  simplified  model  predicts  a  minimum  ignition 
o 

energy  of  ~  8xl02  ergs.  The  full  simulation  does  not  show  ignition, 
but  predicts  that  some  burning  does  occur  and  the  flame  is  eventually 
quenched  (Fig.  9).  Thus  in  the  regime  for  which  both  models  agree,  we 
have  in  fact  tested  them  both.  In  the  regime  where  they  dc  not  agree, 
we  must  then  figure  out  what  physics  is  missing  from  the  similarity 
model.  When  this  is  done,  we  can,  in  effect,  use  the  detailed  model  to 
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Fig.  7  -  Calculated  results  from  the  similarity  solution  plotted  as  a  function  of  time.  Here  Rc  Is  the  char¬ 
acteristic  radius  for  energy  deposition,  A  is  the  nonlinear  amplitude  of  the  temperature  and  density  functions, 
T(R=*0)  is  the  central  temperature  and  I  is  the  induction  parameter.  The  indicates  the  predicted  induction 

time. 


To  =  1.0  x  10  4  sec 
Rn  =  0.1  cm 
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.025  cm 


POSITION  (cm) 

Fig.  9  -  Temperature  as  a  function  of  position  at  four  times  in  the  detailed  flame  model 
The  flame  dies  out  even  though  the  similarity  solution  indicates  it  should  not. 
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build  phenomenology  into  the  similarity  solution.  The  similarity 
solution  has  tested  the  full  detailed  model,  and  the  full  detailed 
model  has  shown  us  the  limits  of  the  similarity  solution. 

V.  Testing  the  Model  with  Experiment 

Although  comparisons  between  analytic  theory  and  model  results  can 
be  used  to  extend  our  understanding  of  the  controlling  processes  in  a 
system  with  limited  physical  complexity,  many  systems  may  preclude  any 
analytic  formulation.  Then  experimental  data  provide  the  only  means  of 
checking  the  accuracy  of  the  model.  Below  we  show  a  case  in  which  the 
results  from  an  experiment  were  used  to  test  a  numerical  model  and  the 
model  results  then  suggested  new  directions  for  the  experiments. 

The  determination  of  the  effects  of  surface  waves  on  submerged 
structures  has  many  practical  applications,  particularly  in  an  ocean 
environment.  Due  to  the  complexity  of  the  problem,  analytic  results  are 
limited  to  idealized  flows  and  geometries.  Amajor  part  of  the  complexity  arises 
from  the  existence  of  the  free  surface  itself.  Not  only  does  the  free 
surface  dominate  the  flow,  but  it  may  become  multiply  connected  when 
sprays,  wave-breaking  or  cavitation  occur.  There  is  usually  no  steady 
state,  and  a  model  for  transient  flow  of  the  fluid  over  the  obstacle 
must  be  used. 

Here  we  describe  the  results  of  a  model  calculation  of  the  wave- 
induced  pressure  forces  on  a  submerged  half -cylinder ,  and  compare  the 
results  with  experimental  data.  The  implications  of  the  comparisons  for 
both  the  validity  of  the  model  and  the  experimental  procedure  will  be 


34 


examined.  Finally,  the  application  of  the  model  to  other  fluid  flows 
and  to  combustion  problems  will  be  discussed. 

Figure  10  illustrates  the  initial  conditions  for  the  numerical 
model.  A  half-cylinder  of  radius  "a"  is  bmerged  in  a  fluid  whose 
undisturbed  free  surface  stands  at  a  height  h  *  2a  over  the  bottom  surface. 
A  progressive  wave  with  wavelength  A  =  5a  is  incident  on  the  cylinder 
from  the  left.  The  sides  of  the  computational  region  are  periodic: 
that  is,  the  physical  system  being  simulated  is  actually  that  of  pro¬ 
gressive  waves  over  a  series  of  half -cylinders.  Periodic  boundary 
conditions  were  chosen  to  avoid  numerical  damping  and  reflection  at  an 
outflow  boundary.  The  calculation  seeks  to  find  the  pressures  at  every 
point  in  the  fluid  as  a  function  of  time,  and  in  particular  the  pressures 
and  pressure  gradient  forces  at  points  on  the  cylinder  surface. 

The  numerical  model  is  based  on  finite  difference  techniques  for 
solving  the  equations  for  inviscid,  incompressible  fluid  flow  using  a 
triangular  grid  which  extends  throughout  the  interior  of  the  fluid  [9^]. 

The  free  surface  and  rigid  boundary  shapes  are  approximated  by  straight 
lines  which  extend  between  points  on  those  surfaces  and  which  define  the 
edges  of  the  computational  grid.  The  governing  equations  are  cast  in  a 
Lagrangian  formalism  so  that  points  originally  lying  on  a  surface  will 
remain  there  at  all  times  during  the  calculation.  Points  interior  to 
the  fluid  will  follow  Lagrangian  pathlines  as  if  they  were  experimental 
marker  particles  in  a  real  fluid.  The  equations  are  differenced  such 
that  vorticity  is  conserved  identically  at  all  times  and  vertex  pres¬ 
sures  are  chosen  to  be  the  values  necessary  to  keep  the  local  fluid 
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volumes  divergence-free.  These  new  pressures  are  in  turn  used  to  advance 
the  velocities  and  update  the  grid  positions. 

The  physical  behavior  of  the  governing  equations  can  be  preserved  in 
the  approximate  difference  equations  being  solved  numerically  by  using  a 
triangular  grid.  A  Lagrangian  grid  will  distort  in  any  non-trivial  flow 
field,  and  as  grid  distortion  becomes  severe  the  calculation  quickly  loses 
accuracy.  However,  a  triangular  grid  can  be  manipulated  locally  in 
several  ways  to  extend  realistic  calculations  of  transient  flows  [9,10]. 
Each  grid  line  represents  a  quadrilateral  diagonal,  and  the  opposite 
diagonal  can  be  chosen  whenever  vertices  move  in  the  flow  to  positions 
which  favor  that  connection.  Such  a  reconnection  involves  just  the  four 
vertices  describing  the  quadrilateral.  No  fluid  moves  relative  to  the 
quadrilateral,  eliminating  one  form  of  numerical  dissipation.  Vertices 
may  also  be  added  or  deleted  to  preserve  the  desired  resolution  by  local 
algorithms  which  involve  only  those  vertices  in  the  vicinity  of  the  grid 
anomaly.  Major  advantages  of  this  technique  are  that  the  algorithms  can 
be  conservative,  they  permit  a  minimum  of  numerical  dissipation  and  yet 
they  require  very  little  computer  time  since  most  of  the  grid  remains 
unaltered. 

Data  for  the  experimental  comparison  was  obtained  through  wave-tank 
experiments  performed  with  a  bottom-mounted  half-cylinder  so  that  pressure 
measurements  could  be  compared  directly  to  the  numerical  results  [11] . 

The  obstacle  was  placed  one-third  of  the  tank  length  from  a  mechanical 
wavemaker  and  at  the  other  end  of  the  channel  a  sloping  porous  beach 
absorbed  95%  or  more  of  the  incident  wave  energy. 
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Results  of  the  experiment  and  of  the  numerical  simulation  are 
shown  in  Figure  11,  together  with  the  results  of  linear  theory.  The 
magnitude  of  the  pressure  fluctuations  as  measured  by  the  experiment  at 
different  points  on  the  cylinder  (9  =  0°  at  the  top  of  the  cylinder) 
are  compared  with  the  predictions  of  the  model.  As  shown  in  Figure  11, 
the  comparison  is  quite  good.  Figure  12  compares  the  calculated  and 
measured  instantaneous  pressure  distribution  around  the  cylinder  for 
the  situation  in  which  the  crest  of  the  progressive  wave  is  near  the 
left  side  of  the  half-cylinder.  Again  the  comparisons  look  good,  but 
now  some  differences  become  evident. 

It  was  found  that  to  within  experimental  error,  all  of  the  observed 
discrepancies  could  be  explained  by  two  factors.  The  first  factor  is 
that  the  model  did  not  exactly  describe  the  physical  situation  in  the 
experiment:  the  wave  tank  had  a  single  cylinder  whereas  the  calculation 
is  for  a  series  of  cylinders.  The  second  factor  was  the  surprising 
result  that  the  roughly  5%  reflected  wave  from  the  wave  tank  signifi¬ 
cantly  affected  the  experimental  results  due  to  modifications  in  the 
dynamic  pressure  fluctuations.  In  this  instance  a  detailed  examination 
of  the  model  and  experimental  results  has  indicated  that  an  experimental 
effect  thought  to  be  small  could  in  fact  cause  noticeable  deviations  in 
the  data  measured. 

The  application  of  this  numerical  technique  to  reactive  flow  is 
relatively  straightforward.  Although  the  example  presented  above  is  for 
homogeneous  flows,  the  extension  to  include  interfaces  involves  no 
basic  changes  to  the  underlying  gridding  scheme,  but  only  the  provision 
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Fig.  11  -  A  comparison  of  numerical  results  with  experiment  and  with  linear  theory  for  the  magnitude 
of  pressure  variations  over  a  half-cylinder  due  to  a  progressive  wave  with  ka  *  w.5,  kh  ■  5. 


pressure  distributions  on  a  half-cylinder  (ka  =  2.5,  kh  =  5),  with  the  wave  crest  just  to 
the  left  of  the  cylinder 


chat  no  interface  sides  are  allowed  to  be  reconnected.  Instead  local 
gridding  anomalies  at  the  interfaces  are  resolved  by  adding  or  subtracting 
vertices  at  Che  interface  and  reconnecting  grid  lines  leading  to  those 
vertices.  That  this  solution  is  viable  is  most  easily  shown  by  Fig.  13 
which  shows  stages  in  the  collapse  of  a  Rayleigh-Taylor  unstable  fluid 
layer  calculated  with  the  same  model.  Here  the  calculation  can  continue 
even  though  the  originally  simply-connected  lighter  fluid  performs  a 
transition  to  a  multiply  connected  fluid  which  includes  "bubbles"  which 
have  been  entrained  by  the  heavier  fluid.  Of  course,  for  reactive  flow 
calculations  a  new  model  would  have  to  be  constructed  based  on  these 
techniques  which  used  instead  the  equations  governing  compres¬ 

sible  fluids  and  which  contained  the  added  chemical  and  thermal  equations 
besides . 

VI .  Conclusion 

Deatiled  modelling  of  laminar  reactive  flows,  even  in  fairly  compli¬ 
cated  geometries,  is  certainly  well  within  our  current  capabilities.  In 
this  paper  we  have  shown  several  ways  in  which  these  techniques  may  be 
used.  As  the  physical  complexity  we  wish  to  model  increases,  our 
footing  becomes  less  sure  and  more  phenomenology  must  be  added.  For 
example,  we  might  have  to  add  evaporation  laws  at  liquid-gas  interfaces  or 
less  well-known  chemical  reaction  rates  in  complex  hydrocarbon  fuels. 

Perhaps  the  biggest  problem  facing  combustion  modelling  now  is 
turbulence:  there  are  no  excellent  or  even  good  methods  of  including  such 
effects  in  our  calculations.  At  best  we  have  a  number  of  phenomenological 
models  with  limited  ranges  of  validity  and  which  imply  a  steady  state. 
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We  believe  that  devising  a  way  to  handle  this  difficult  problem  of 
strongly  coupled  multiple  time  and  space  scales  is  the  challenge  we 
currently  face. 
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